Metastable liquid-liquid phase transition in a single-component system with only one 

crystal phase and no density anomaly 



(N 
O 
O 
(N 

D 
00 



o 



o 

m 
> 

m 
(N 



I 

O 

o 



X 



G. Franzese^'^j] G. Malescio^, A. Skibinsky\ S. V. Buldyrevi, and H. E. Stanley^ 
^ Center for Polymer Studies and Department of Physics, Boston University, Boston, MA 02215, USA 
^ Dipartimento di Ingegneria dell'Informazione, Seconda Universita di Napoli, 
Istituto Nazionale Fisica della Materia UdR Napoli and CG SUN, 1-81031, Aversa, Italy 
Dipartimento di Fisica, Universita di Messina and Istituto Nazionale Fisica della Materia, 1-98166 Messina, Italy 

(Dated: February 1, 2008) 

We investigate the phase behavior of a single-component system in 3 dimensions with spherically- 
symmetric, pairwise-additive, soft-core interactions with an attractive well at a long distance, a 
repulsive soft-core shoulder at an intermediate distance, and a hard-core repulsion at a short dis- 
tance, similar to potentials used to describe liquid systems such as colloids, protein solutions, or 
liquid metals. We showed [Nature 409, 692 (2001)] that, even with no evidences of the density 
anomaly, the phase diagram has two first-order fiuid-fiuid phase transitions, one ending in a gas- 
low-density liquid (LDL) critical point, and the other in a gas-high-density liquid (HDL) critical 
point, with a LDL-HDL phase transition at low temperatures. Here we use integral equation cal- 
culations to explore the 3-parameter space of the soft-core potential and we perform molecular 
dynamics simulations in the interesting region of parameters. For the equilibrium phase diagram we 
analyze the structure of the crystal phase and find that, within the considered range of densities, 
the structure is independent of the density. Then, we analyze in detail the fiuid metastable phases 
and, by explicit thermodynamic calculation in the supercooled phase, we show the absence of the 
density anomaly. We suggest that this absence is related to the presence of only one stable crystal 
structure. 

PACS numbers: 61.20.Gy, 65.20.+W, 64.70.Ja, 64.60. My 



I. INTRODUCTION 



Soft-core potentials have been widely used to study a 
variety of systems such as liquid metals, metallic mix- 
tures, electrolytes and colloids, as well as anomalous liq- 
uids, like water and silica |], |, |, |, |, |, |, |, ||, |l^, 
|T2| , |l3| , p^ , p^ . In these models, the specific structural 
characteristic at the molecular (or atomic) level are ne- 
glected and the molecules (or atoms) are represented by 
simple spheres. Quantum effects (such as the quantum 
nature of chemical interactions) and classical effects (such 
as the Coulomb interaction) are modeled through a phe- 
nomenological isotropic pair potential. The advantages 
of this approach are that while these potentials are simple 
enough to be treated analytically they still allow a 
qualitative comparison with the experiments. Moreover, 
they can be studied by means of numerical simulations 
less time-consuming than those of realistic models . 

We consider an off-lattice model in three dimensions 
(3D) jll] related to the soft-core potentials studied by 
Hemmer and Stell |^ for solid-solid critical points. Our 
model shows a phase diagram with two fluid-fluid phase 
transitions, a feature recently seen in experiments on 

t)sphorus and conflrmed by speciflc simulations 
. In Ref. ||11|| we showed that both first-order fluid- 
fluid phase transitions end in critical points, a low-density 
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critical point Ci, and a high-density critical point C2. 
For the considered potential, both transitions occur in 
the supercooled phase with respect to the crystal phase. 

The hyp othesis of a second critical point has been pro- 
posed po| as a way to rationalize the density anomaly — 
i.e., the expansion upon isobaric cooling — in network- 
forming fluids, such as water ||l|, |2^, |2^, carbon [ p2[ , 
silica and silicon 12^. Consequently, both the 
experimental l 25i |2^, |24" ^^, p9[ and the theoretical 
@, § |, |l|r|l3~|l|, jsi) investigations about the 
possibility of a second critical point have been focused 
on systems with the density anomaly (Fig. ^). How- 
ever, the results of Ref. ||ll[ have shown that the presence 
of the critical point C2 does not necessarily induce the 
density anomaly, indicating that the quest for simple liq- 
uids with two critical points is not restricted to systems 
with densities exhibiting anomalous behavior (Fig. |l|b). 
In this paper we push forward the analysis, by studying 
the equilibrium phase diagram and showing that the sys- 
tem introduced in Ref. has one single crystal phase, 
in the range of considered densities, suggesting that the 
absence of density anomalies is related to the presence of 
only one stable crystal structure. 

To reach this goal we organize the work in the fol- 
lowing way. After the definition of our soft-core poten- 
tials (Section II), (i) we describe in detail the integral 
equations in the hyper-netted chain (HNC) approxima- 
tion (Section III. A) and (ii) we generate new calculations 
using different assigned parameters for the pair potential, 
thus establishing a bridge between the potential studied 
in Ref. jll] and the potential investigated in Ref. [ 10[ 
and rationalizing how — in addition to the critical point 
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FIG. 1: Schematic pressure-temperature P — T phase dia- 
grams with two critical points. Solid lines represent first or- 
der phase transition lines, circles represents critical points, (a) 
Phase diagram with critical point C between the gas and the 
uniform liquid phase at high T and low P and critical point 
C' between the low-density liquid (LDL) and the high-density 
liquid (HDL) at low T and high P. This phase diagram has 
been proposed for water and has been shown to be consistent 
with the density anomaly, (b) Phase diagram with critical 
point Ci between the gas and the LDL at low T and low 
P and the critical point C2 between the gas and the HDL 
at high T and high P. Increasing P at constant tempera- 
ture Ta below Ci (dashed line at To), the system undergoes 
a first order phase transition between the gas and the LDL 
phase, followed by a first order phase transition between the 
LDL and the HDL phases. Increasing the pressure at con- 
stant temperature Tb above Ci but below C2 (dashed line at 
Ti,), the system undergoes only a first order phase transition 
between the gas and the HDL. The square represents the gas- 
LDL-HDL triple point. This phase diagram has been found 
in Ref. with no evidence of the density anomaly. 



Ci — the critical point C2 arises as a function of the pa- 
rameters of the pair potential (Section III.B). Then, (iii) 
we describe in detail the molecular dynamics (MD) simu- 
lations, studying the equilibrium phase diagram and find- 
ing that the only stable crystal structure, in the range 
of simulated densities, has two characteristic lattice dis- 
tances (Section IV. A), (iv) By MD simulations we ana- 
lyze also the supercooled liquid phase and the metastable 
fluid- fluid phase transition (Section IV. B); (v) in order 
to construct an accurate phase diagram and to study the 
finite size effect, we perform new MD simulations in ad- 

(Section 
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dition to the calculations presented in Ref. 
IV. C). Hence, (vi) we study the radial distribution func- 
tion, comparing the HNC predictions with the MD re- 
sults and analyzing the composition of the system within 
the fluid- fluid coexisting regions (Section V). Finally, 
(vii) we address the density anomaly issue by present- 
ing the explicit thermodynamic calculation, based on the 
MD phase diagram, that allows us to exclude the pres- 
ence of the density anomaly (Section VI). We give our 



conclusions in Section VII. 



II. SOFT-CORE POTENTIALS 

Among the isotropic potentials, much attention has 
been devoted to soft-core potentials, which have a finite 
(soft-core) repulsion added to the infinite (hard-core) re- 
pulsion. The infinite repulsion is due to the impenetra- 
bility of the spheres. The finite repulsion represents the 
combination of all the quantum and classical repulsive ef- 
fects averaged over the angular part. It has been shown 
Q that a weak effective repulsion can be derived by a 
first principle calculation for liquid metals js). To un- 
derstand the possibility of the solid-solid critical point in 
such material as Ce and Cs, Hemmer and Stell [|| pro- 
posed a soft-core potential with an attractive interaction 
at large distances, performing an exact analysis in ID. 
Over the past thirty years, several other soft-core poten- 
tials with one or more attractive wells were proposed and 
studied with approximate methods, or numerical simula- 
tions in 2D, to rationalize the properties of liquid met- 
als, alloys, electrolytes, colloids and the water anomalies 
0, |, ||, 11 J I |ll|, [||l|, |14[ |l| . It has been re- 
cently shown [|nj, for the first time in 3D with MD, that 
an appropriate soft-core potential with an attractive well 
is able to show two supercooled liquids of different den- 
sities, with two critical points. Similar results have been 
reported for a soft-core potentials with a linear repulsive 
ramp by Monte Carlo simulations . 

Following Ref. , we define the isotropic pair poten- 
tial U{r), as a function of the pair distance r (inset in 
Fig.|): 



C/(r) = 




for r < a 

for a < r < b 

for b < r < c 

for c < r , 



(1) 



where a is the hard-core distance and b is the soft-core 
distance. For a < r < 6, the spheres interact with a 
finite (soft-core) repulsive energy Ur. For b < r < c, 
the pair's interaction is attractive with energy —Ua < 
0. The distance c is the cut-off radius beyond which 
the pair's interaction is considered negligible. For sake 
of comparison with Ref. [|lo|, |l5|, we will consider both 
Ur>0 and Ur < 0. 

The step-like shape in Eq. (|l]) has the advantage of 
being defined by only three parameters: the width of 
the soft-core in units of the hard-core distance wrJ a = 
[b — a) /a, the width of the attractive well in the same 
units wa/cl = {c—b)/a and the soft-core energy in units of 
the attractive well energy Ur/Ua- To explore the phase 
diagram of the model as a function of these three param- 
eters in an approximate, yet fast way, we use the integral 
equations in the HNC approximation. 
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III. THE INTEGRAL EQUATIONS IN THE 
HYPER-NETTED CHAIN APPROXIMATION 

In this section we present details of the integral equa- 
tions and the HNC approximation adopted in Ref. Q 
and in this work. The radial distribution function g(r) 
plays a central role in the physics of fluids This 
quantity is proportional to the probability of finding a 
particle at a distance r from a reference one and is the 
ratio of local to bulk density at distance r, with 
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n(r) 



(2) 



where n(r) is the number of particles at a distance be- 
tween r and r + dr from the reference particle and p is the 
number density, assumed to be independent of r (uniform 
system). The radial distribution function goes to 1 for 
large r and is always 1 for a random spatial distribution 
of particles. To represent deviations from randomness, 
the total pair correlation function h(r) = g{r) — 1 is in- 
troduced. 

These functions are relevant because they are directly 
measurable by radiation scattering experiments and are 
related to the thermodynamic properties of the fluid. A 
fundamental relation between structure and thermody- 
namics is given by 



ksTpKx = 1 + p h(f) dr 



(3) 



where = {dp/dP)T/ p is the isothermal compressibil- 
ity, P is the pressure, and fc^ is the Boltzmann constant. 
Provided that the particles interact through pairwise- 
additive forces, other thermodynamic properties of the 
fluid — such as the internal energy — can be expressed us- 
ing integrals over the pair correlation function. 

The function h{r) is the result of the interaction of 
all the particles in the system. Formally, h{r) can be 
decomposed into (i) the contribution coming from the 
direct interaction between two particles at distance r, 
called c(r), and (ii) the contribution due to the indirect 
interaction propagated through any other particle in the 
system. This second contribution is written in turn as an 
integral convolution of direct correlations and total pair 
correlation. 

This decomposition, for uniform systems, is expressed 
by the Ornstein-Zernike relation 



h{r) = c(r) + p c{r')h{\r — r'\) dr' 



(4) 



Equation (H) is also the formal definition of c(r). Both 
h(r) and c(r) in Eq. are unknown functions, thus to 
solve this equation, one needs another relation (closure) 
between these two functions. This relation is provided 
by the diagrammatic expansion of g{r) which, after 
formal summation, yields the functional relation 

gir) = exp[-/3t/(r) + h{r) - c{r) + d{r)] , (5) 



where U{r) is the inter-particle potential, (3 = l/^ksT) 
and d{r) is the sum over a specific class of diagrams 
[bridge diagrams) [ p^ . Since d{r) cannot be calculated 
exactly, one resorts to approximate expressions. The sim- 
plest approximation assumes d{r) — (HNC closure) 
|32| . One expects this approximation to work better 
at lower p, where the direct correlation function c(r) is 
more relevant then the correlation propagated through 
the other particles. However, our results (see Section V) 
will show that this intuitive observation is not straight- 
forward, at least for soft-core potentials. 



A. The iterative procedure 



The solution of the integral Eqs. (^ with the HNC 
closure is obtained through a numerical iterative proce- 
dure whose essential scheme is the following. Under the 
assumption d{r) = 0, one can write g{r) as 



g{r) = exp[-/3C/(r) -I- 0{r)] 



(6) 



where the function 9{r) = h{r) — c(r) has the remarkable 
property of being a continuous function of r, even for 
discontinuous potentials (as in this paper). From the 
definitions of h{r), 9{r) and Eq. (||), one can derive the 
equation 



c(r) = exp[-/3C/(r) + e{r)] - 0{r) - 1 



(7) 



By using the Fourier transform f(q) = / f{r) exp(ig-r) df 
defined for a generic function /{f), from Eq. (|[) we obtain 



h{q) = c{q) + p c{q) h{q) . 
Or, using the definition of 0{r), we have 



Oiq) 



cHq) 



1 - p c{q) 



(8) 



(9) 



The numerical iteration is based on Eqs. (0) and (^. 

We start by choosing an initial guess for 9{r). A rea- 
sonable input, at least at high temperatures, is the 9{r) 
of a fluid of hard spheres with diameter a. In fact, at high 
temperatures our potential can be approximated with a 
simple hard-core repulsion. We can calculate the corre- 
sponding 0{r) by making use of the Percus-Yevick inte- 
gral equation [|6| which for hard spheres can be solved 
analytically. Next, at constant p, we decrease the tem- 
perature of ST and we perform calculations at fixed p 
and T by using as input the 9{r) obtained as solution at 
p and T + ST. 

From the chosen guess of 9{r) we calculate c(r) by us- 
ing Eq. (^. Its Fourier transform c{q) is used in Eq. (||) 
to calculate 9{q). Its inverse Fourier transform provides 
a new 9{r) that is used as a new input for the next cycle. 
We evaluate the functions on M = 2048 discrete points 
r„ = mSr^ with m = 1, . . . , M and Sr = O.Ola. Succes- 
sive iterations of the elementary cycle defines a succession 
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B^'') (r) , where fc = 1, 2, . . . is the number of the iteration. 
If the difference between two consecutive elements of this 
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succession, 



1/2 



(10) 



decreases for increasing fc, the succession converges to- 
wards a 0*{r) that is solution of our integral equations. 
The iteration process is stopped when A < 10"''. 

Based on this iterative procedure, different algorithms 
can be used to improve the accuracy and rapidity of con- 
vergence of the numerical solution of HNC equations. 
However, independent of the algorithm used, there ex- 
ists a region in the p-T plane where no solution can be 
found, i.e., for any p, there is a T below which the numer- 
ical algorithm does not converge, defining an instability 
line in the p-T plane. 



B. The HNC instability line 

The nature of the locus of instabilities of the HNC 
equation and its relationship with the spinodal line of 
the fluid was investigated for a hard-core potential plu s 
an attractive Yukawa tail in a number of papers |4| . 
These studies showed that the isothermal compressibility 
does not diverge as the temperature is lowered and the 
instability region is approached from above. This conclu- 
sion was definitively assessed through extensive numeri- 
cal calculations both for the hard-core Yukawa fluid 
and other model potentials, showing that this behavior 
is directly correlated to the existence of multiple HNC 
solutions. The analysis developed was based on a careful 
treatment of the low-fc behavior of the Fourier transforms 
of the correlation functions required by the iterative pro- 
cedure. A further theoretical support to these results was 
given by an analysis on models for an ionic fluid and 
a monoatomic Lennard-Jones fluid. 

In the light of the above mentioned studies, an identifi- 
cation of the instability line of the HNC equation with the 
spinodal line of the fluid, which is characterized by a di- 
verging compressibility, is not possible. Keeping in mind 
this limitation, one can nevertheless observe that for a 
large number of simple fluid pair potentials the shape of 
the instability line qualitatively resembles the region of 
spinodal decomposition of the fluid. Also for our poten- 
tial the comparison of the HNC calculations with the MD 
results (Section ^ shows that the HNC instability line is 
qualitatively consistent with the spinodal line |^^. Thus, 
studying the modifications of the instability line as the 
potential parameters are varied can yield some approxi- 
mate, yet useful, informations on the phase behavior of 
the fluid. 
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FIG. 2: Instability line of the HNC equations in 3D in the 
p-T plane for the pair potential in Eq. (|l|), with wr/o = 0.4, 
wa/o, = 0.3 and (from top to bottom) Ur/Ua = —0.5, 0.0, 
0.5. For Ur/Ua = —0.5, the pair potential recovers the one 
studied in ID and 2D in Ref. [ 10|, [l^] . The symbols represent 
the calculations and the lines are guides for the eyes. Inset: 
the isotropic pairwise-additive potential U(r-) as a function of 
the inter-particle distance r; a is the hard-core distance, b is 
the soft-core distance, c is the maximum attractive distance, 
— Ua < is the attractive energy and Ur is the repulsive 
energy. 



C. The results 

First, we calculate the instability line of the HNC equa- 
tions for the potential investigated in Refs. ||l^, |l^. The 
corresponding parameters are wa/a = 0.4, wa/u — 0.3 
and Ub./Ua = —0.5. In this case, the soft-core is given by 
two attractive wells with different depths. Calculations 
in ID and 2D |l^, |l^ have shown a water-like density 
anomaly. Therefore, it is interesting to analyze the phase 
diagram in 3D. However, the instability line for this case 
(Fig. ^ is similar to the spinodal line usually exhibited 
by a simple fluid, e.g. interacting via a Lennard-Jones 
potential with the maximum of the spinodal line corre- 
sponding to the liquid-gas critical point. Upon increasing 
Ur/Ua to 0.5 (Fig. ||), the only evident change of the in- 
stability line is a shift toward a lower T as a result of the 
overall decrease of the inter-particle attraction, with no 
hints of a second critical point. A small shift to lower p 
is also seen. This behavior is more evident for larger wji. 

Next, we consider a potential with a larger wr 
(wji/a = 0.7, wa/cl = 0.3). The instability line is calcu- 
lated for several values of Ur/Ua (Fig- H)- Upon increas- 
ing Ur/Ua, we now find not only the shift to a lower 
T, but also an evident shift of the maximum of the line 
(i.e. the critical point, assuming that the instability line 
represents the behavior of the spinodal line) to a lower p. 
This result can be rationalized by observing that, passing 
from Ur < to Ur > 0, the soft-core becomes more and 
more difficult to penetrate and the system passes from a 
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WA/a = 0.2. The results (Fig. |) show that for 0.4 < 
Ur/Ua < 0.6, the instabihty hne has two well-distinct 
local maxima, suggesting the possibility of two critical 
points in the phase diagram for the fluid phases . For 
Ur/Ua < 0.3 or Ur/Ua > 0.7, the instabihty hne shows 
just one maximum, similar to the typical spinodal line of 
a fluid of hard spheres with diameter a or 5, respectively, 
attracting via a square well of width wa- As a conse- 
quence of this analysis, we choose wr/u = 1, wa/cl = 0.2 
and Ur/Ua — 0.5 as the set of parameters for the poten- 
tial used in the MD calculations in 3D llTll. 



IV. THE MOLECULAR DYNAMICS 
APPROACH 



FIG. 3: Instability lines, as in Fig. ^ for the pair potentials 
with parameters wr/u = 0.7, wa/o, = 0.3 and (from top to 
bottom) Ur/Ua = -0.5, 0.0, 0.5, 0.7, 1.0. 
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FIG. 4: Instability lines, as in Fig. ^ for the pair potentials 
with parameters wr/u — 1.0, wa/cl = 0.2 and (from top to 
bottom) Ur/Ua = 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7. The 
full symbols correspond to the set of parameters selected for 
the MD calculations. 



potential with a hard-core a and an effective attractive 
range wa + wr to a potential with an effective hard-core 
6, for Ur/Ua large enough, and an attractive range wa- 
As a consequence of the increase of the effective hard- 
core, the critical density decreases and, as a consequence 
of the decrease of the effective overall attraction, the crit- 
ical temperature decreases. 

Comparing Fig. |^ and Fig. |^, we notice an important 
difference. In the case of larger wr, (Fig. as Ur/Ua 
increases, the temperature of the instability line does not 
decrease for increasing p, but becomes rather flat. This 
result suggests that the instability line might develop a 
second maximum at larger values of p for even larger wr . 

We thus consider a potential with wr/g = 1.0 and 



In this Section we give extensive details on the MD 
method and we extend the analysis performed in Ref . 11 1, 
including new calculations for the crystal phase, the crys- 
tal nucleation process and the metastable phases. We 
perform MD simulations at constant number of particles 
N of unit mass m, at constant volume V, with periodic 
boundary conditions and at constant average tempera- 
ture T. We present the results for N = 490, 720 and 
N — 1728. The average temperature is set by coupling 
the system to a thermal bath at the assigned T, with 
a thermal exchange coefficient per particle between the 
system and the bath equal to A: = 0.015 {UA/mY^^kB/a. 
We use a standard collision event list algorithm j3^] to 
evolve the system and a modified Berendsen method to 
achieve the desired T |Q . 

The pressure is calculated by using the virial expres- 
sion for a step potential ll^ 



/ ^ 1 

^ TO / >. - , 1 



N 



with ^ - J sum over the particles pairs undergoing 
a collision in the time interval At = (lO^ma^ /Ua)^^'^ , 
hereafter used as unit of time, and with Avi = v'^ — Vi, 
where Vi and v'^ are the velocities of the particle i at 
position fi before and after the collision with particle j 
at position rj. 



At ^ 



'Avi ■ ifi - rj) 



(11) 



A. The crystal 

First, to locate the equilibrium crystal line we simulate 
a crystal seed surrounded by the gas. We prepare a crys- 
tal seed by cooling at T = 0.45C/yi/fcs a gas configuration 
with density p = 0.018/a^. 

The crystal (Fig. |^) is the effect of the competition 
between the hard-core repulsion at distance r = a and 
the attraction at distance r — b. The resulting struc- 
ture is reminiscent of the close packing of hard spheres 
with diameter a or b, but the competition gives rise to 
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FIG. 5: The MD configuration equilibrated at ksT/UA = 
0.45 and a^p = 0.018. Darker particles are farther away 
from the observation point, (a) The crystal, with defects, 
surrounded by the gas. Bonds connect particles at distance 
1 < r/a < 1.2. The radius of the particles is not in scale with 
the distances. A typical ring of 8 particles (octagon) is plot- 
ted with a larger radius (b) A section of the crystal. Bonds 
connect particles at distance 2 < r/a < 2.2. The radius of 
the particles is in scale with the distances. Note, in the up- 
per part of the panel, a ring of 12 nn particles (dodecagon), 
connected by bonds to two central particles, (c) The same 12 
particles of the section in (b) are plotted with a larger radius 
with respect to the other particles. Bonds are like in panel 
(a), (d) A rotation of 40 degrees around a central horizontal 
axis of the section in (c) reveals the 8-fold symmetry observed 
in (a). 



new symmetries (Fig. The minimum in the inter- 
particle interaction potential at b < r < c would induce 
a face-centered-cubic crystal with lattice space ranging 
from & to c and a characteristic six-fold symmetry on one 
projection plane with seven particles to form a triangu- 
lar lattice. However, due to the soft-core, the system 
can allocate particles at the hard-core distance a = 6/2. 
This induces a twelve-fold symmetry, placing an aver- 
age number of twelve, almost on-plane, nearest neighbor 
(nn) particles at a distance 1 < r/a < 1.2 to form a 
dodecagon around two particles. These two nn parti- 
cles are next nearest neighbors to the dodecagon, at a 
distance 2 < r/a < 2.2, and are placed on a line al- 
most perpendicular to the plane individuated by the do- 
decagon (Fig. ||b). This structure is distorted in such a 
way to form non-closed chains of nn particles, that wrap 
along another axis to give rise to an eight-fold symmetry 
(Fig.|c,d). 

By analyzing the crystal structure obtained from the 
MD simulations, we conclude that the position of the 
particles in the crystal can be described by r = i ■ a + 
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FIG. 6: The artificial crystal configuration equilibrated for a 
time 10^ At with a MD simulation at T = 0.03i7A/fcs. Bonds 
connect particles at distance r/a < 1.2. The radius of the par- 
ticles is not in scale with the distances. Greater particles are 
closer to the observation point. The configuration contains 
15 cells. The central cell is emphasized by darker bonds, (a) 
Each cell has 4 particles at the corners of a rectangle (rvi 
with m = 1, ... 4 in Table whose long sides form two tri- 
angles with two particles on the same plane {f,n with m = 5, 6 
in Table and the short sides forms two tetrahedra, each 
with two more particles (fvi with m = 7, ... 10 in Table ^ . 
(b) The crystal configuration rotated by 7r/4 around a central 
horizontal axis shows the 8-fold symmetry seen in Fig. |^ and 
d. (c) A further rotation of 7r/4 around the same axis shows 
the dodecagons seen in Fig. plb and c. (d) A rotation of it/2 
around a central vertical axis shows again the octagons seen 
in Fig. Ha and d. 



TABLE I: The coordinates of the lattice vectors as obtained 
after an equilibration time 10^ At at T = 0.03C/a/A;b of the ar- 
tificial crystal (Fig. ^ proposed to describe the crystal struc- 
ture found in the MD simulations (Fig. The errors on the 
parameters are on the last decimal digit and decrease as the 
square root of the time of averaging. 
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a/a 


1.95 


0.00 


0.11 


b/a 


O.OO 


3.41 


0.00 


c/a 


O.OO 


0.00 


1.94 



k ■ c + fmi where r^, for m = 1, . . . 10, are the 



coordinates, with respect to the center of the cell, of the 
10 particles forming a crystal cell, a, b and c are the 
lattice vectors describing the position of the center of the 
cell and i, j and k are integers such that i + j + k is 
even. We estimate the lattice vectors (Table ^ and the 
coordinates of the particles forming the cell (Table ||) 
after an equilibration time 10^ At at T = 0.03UA/kB for 
of an artificial crystal placed in the vacuum (Fig. o). The 
resulting density of the crystal is a^p ~ 0.39. 

Surface effects could be responsible for the tilt that can 
be seen in Fig. Bd. In a system with N = 720, this tilt 



TABLE II: The coordinates 



for 



1, ... 10 of the 10 particles forming a crystal cell, with respect 
to the center of the cell, as obtained after an equilibration 
time 10^ At at T = 0.03UA/kB of the artificial crystal (Fig.||). 
The characteristic distances, with an error on the last decimal 
digit, are: h/a = 0.53, h/a = 0.59, h/a = 0.07, k/a = 1.43, 
k/a = 0.05, hja = 0.03, hja = 1.40, hja = 0.52. The errors 
decrease as the square root of the time of averaging. For each 
particle m of the cell, we denote with the number of par- 
ticles in the crystal at a distance r < b (in the soft-core) and 
with the number of particles at a distance b < r < c (in 
the attractive well). 
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FIG. 7: The artificial crystal configuration for A'^ = 720 par- 
ticles equilibrated with a MD simulation at T = 0.10(7A/fcfl 
(a) and at T = 0.52(7^/^3 (b). Bonds connect particles at 
distance r/a < 1.2. The radius of the particles is not in scale 
with the distances. The crystal seeds are in equilibrium with 
the gas phase and show many defects. The tilt present at the 
lower T in (a), disappears at the higher T in (b). 



disappears when the sample is equilibrated at higher T 
(Fig. 0). 

We compare the g{r) (Fig. |^) of the MD crystal in 
Fig. H and of the artificial crystal in Fig. both equili- 
brated at T = OASUA/kB- Both functions show peaks 
located at the same distances, with two large peaks at 
r/a = I and r/a — 2, consistent with the presence of the 
two characteristic distances a and b in the potential. The 
comparison confirms that the proposed crystal is a good 
representation of the crystal structure generated by the 
MD simulation. The slightly different intensities of the 
peaks of the g{r) of the two systems are probably due to 
the defects of the MD crystal. 

The validity of the artificial crystal as a good descrip- 
tion of the real crystal structure is confirmed also by 
the evolution of the potential energy per particle (in- 
set in Fig. 0) when the MD crystal and the artificial 
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FIG. 8: The radial distribution functions g{r) for the MD 
crystal (solid line) and the artificial crystal (dotted line), both 
equilibrated at T = OASUA/ks , are very close to each other. 
Inset: the potential energy density U/N for the MD crystal 
(solid line) with A'^ = 490 and the artificial crystal (dotted 
line) with A" = 720, both equilibrated at T = 0.6017^/^3 
starting from a configuration equilibrated at T = 0.48(7^/^3, 
as function of the time divided by the number of particles A'. 



crystal are heated, from the configuration equilibrated 
at T = 0A8UA/kB, to T = 0.60UA/kB. Both samples 
equilibrate to the same energy. The starting potential 
energy is, as expected, in both cases greater than the 
ground state energy Uq/N = —8A5Ua, calculated from 
the number of particles at distance a < r < b and at dis- 
tance b < r < c (Table ||), due to surface effects. We find 
analogous results for the evolution of the kinetic energy. 

To test if the system has more than one crystal struc- 
ture as function of the density, we cool dXT — Q.6UA/kB 
a fluid configuration equilibrated at T = 0.8UA/kB and 
p = 0.267/a'^, and compare the resulting g{r) with 
the case of the crystal seeds at T = 0A5UA/kB and 
p = 0.018/a'^, finding no relevant differences (Fig. ^). At 
the same time, the attempt of finding alternative artifi- 
cial crystal structures has revealed, after an appropriate 
equilibration, that the only stable structure is the one 
presented in Fig. ^. We therefore assume that the sys- 
tem, at least for this choice of the potential's parameters, 
has one single crystal structure independent of p, within 
the considered range of densities. 

Starting from a configuration with the crystal seed 
described above, we equilibrate the system at different 
densities and temperatures. We define the system to be 
in the solid phase if, after a time 10^ {ma'^ /UaY^'^ — 
3 X 10'^ Ai, the crystal seed is growing, or we consider 
it in a fluid phase if the seed is melting. The cases in 
which the trend is not clear within the simulation time 
are considered as belonging to the first-order transition 
region pl| . The crystallization pressure rapidly increases 
with p and T, giving a first-order transition line (in the 
thermodynamic limit) that separates the equilibrium P- 
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FIG. 9: Comparison of the radial distribution function g{r) 
calculated for two MD configurations obtained by cooling the 
system at different densities: the solid line is for the config- 
uration at ksT/UA = 0.60 and a'^p — 0.267, the dashed line 
for the configuration at kBT/UA = 0.45 and aP' p — O.OfS. 
The two functions are very close to each other for distances 
r/a < 4, showing that the crystal structure is the same at 
both densities. The difference between the two functions is 
consistent with the presence of defects and of the surrounding 
gas. 



p phase diagram in a high-T fluid and a low-T crystal 
(Fig. 0) p. 



B. The supercooled liquids 

At equilibrium, there is no (stable) liquid phase. A 
phase diagram without equilibrium liquid phase is ex- 
pected for systems with an inter-particles poten- 
tial with a narrow attractive part, such as the one we 
are considering here. However, the liquid is present as a 
metastable (supercooled) phase with respect to the crys- 
tal phase Q. To study the metastable phase diagram, 
we equilibrate the system for each p from a gas config- 
uration at T = 0.70UA/kB and then rapidly cool it to 
the desired T > 0.57J7^/fc_B, calculating P, g{r) and the 

total potential energy U = J2i<j ~ '^iD- 

We find that the supercooled fluid phase has a life- 
time longer than 3 x 10^^At (the standard length of our 
simulations) for p < 0.20/a^ at T sa 0.57UA/kB, for 
p < 0.27/a^ at T « OMUA/ks and for p < OM/a^ at 
T « Q.70UA/kB- The system is equilibrated in the fluid 
phase for t > 20At, after which we average P, g{r) and 
U over the time. We calculate each state point by aver- 
aging the conflgurations for 3 x 10^ At < t < 3 x l{filS.t. 
We estimate the errors by dividing the configurations in 
90 non-overlapping intervals of 30At, that we assume in- 
dependent. 

For larger p, the system spontaneously crystallizes (ho- 
mogeneous nucleation process). Thus, we only average 
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FIG. 10: The MD P-p phase diagram. The thick dashed 
line is the gas-crystal first-order transition line at the equi- 
librium. The calculations in the region below this line are 
for the metastable fiuid states. Main panel: the diamonds 
(full and open) are the MD calculations for (bottom to top) 
ksT/UA = 0.570, 0.580, 0.590, 0.600, 0.610, 0.620, 0.630, 
0.640, 0.650, 0.660, 0.670, 0.675, 0.685, 0.700. The dotted 
lines are the isotherms calculated by polynomial interpola- 
tions of the points at constant T and, at the same time, of 
the points at constant p. The circle at p — 0.1/ is the 
gas-LDL critical point Ci. The square at p = 0.306/a'^ is 
the gas-HDL critical point C2. The solid thick lines connect- 
ing the local minima and maxima along the isotherms are the 
spinodal lines associated to each critical point and the shaded 
regions are the associated mechanically unstable regions. The 
dashed lines, passing through the critical points, are the coex- 
istence regions associated to each critical point. The meeting 
point of the gas-HDL coexistence line with the gas-LDL co- 
existence line gives the possible triple point (full triangle at 
p ~ 0.12/a^). Where not shown, the errors are smaller then 
the symbols. Inset: enlarged view around Ci. Symbols are as 
in the main panel. The diamonds are the MD calculations for 
(bottom to top) kgT/UA = 0.580 0.590, 0.595, 0.600, 0.603, 
0.606, 0.609, 0.620, 0.630, 0.640, 0.660. 



over configurations that occur before nucleation. To be 
certain that our estimates are carried out in the fiuid 
phase, we study 



S{q,t) 




Jq-lrj{t)-fk(t)] 



(12) 



where fj (t) is the position of particle j at time t and q 
is the wavevector. At equilibrium, the average of S{q,t), 
over the time and the wavevectors with the same module, 
is the structure factor <5'(g), describing the spatial correla- 
tion in the system. Therefore, S{q,t) describes the time- 
evolution of the spatial correlation along the wavevector 
q. In particular, for a crystal-like configuration, with a 
long-range order, there is at least one wavevector such 
that S{q,t) 0{N), while for a fiuid-like configuration 
is S{q,t) - 0(1) for afl q. 

The time-evolution of S{q., t) for a typical simula- 
tion inside the nucleation region is presented in Fig. |l^. 
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FIG. 11: MD calculations, (a) Time evolution of the struc- 
ture factor S{q,t) at T = 0.62UA/kB and p = 0.27/a^, for 
wavevectors with modulus q — 1/a and for a time t/At — 
3000; in the time interval "A" with 200Af < t < 500At, it 
is S{q,t) ^ 0(1) for any q; in the time interval "B" with 
700At < t < ISOOAt, for six wavevectors there is an in- 
crease in S(q = l/a,t); in the interval "C" with 2100Af < 
t < SOOOAt, for the same six wavevectors there is a larger 
increase, (b) As in (a) but for g ~ 12/a ~ 47r/a; in this case 
there is a large increase in S{q ~ 12/a,t), more then one or- 
der of magnitude, only in the time interval "C" for several 
wavevectors, revealing the formation of a crystal seed, (c) 
The structure factor S{q), given by the average over the di- 
mensionless wavevectors aq with the same modulus and the 
average over the time intervals "A", "B" and "C" of S{q,t). 
The curves for "B" and "C" are offset by 1.5 and 3, respec- 
tively. All the curves go to 1 for large q. In the interval "A" , 
S{q) is liquid-like. In the interval "B", S{q) is still liquid-like 
but with an increase for q ~> 0, while in the interval "C" is 
solid-like, with two large peaks at g ~ 27r/(6/2) — 27r/a and 
q ~ 27r/(a/2), corresponding to the soft-core radius and the 
hard-core radius, respectively, and a large value for q ~* 0. 



FIG. 12: Inset: projections of the 3D configuration of A'^ = 
490 particles in a system of size yi/^ = 12.25a and cor- 
responding to the largest peak in the time interval "B" in 
Fig. [l]]a. The projections are (from left to right) Z vs Y, Z 
vs X and X vs Y. The histograms of number of particles as 
functions of the abscissa is over-imposed on each projection. 
Projections and histograms are shifted for clarity. Each his- 
togram bin corresponds to half of the box size. The dashed 
line shows the average number of particles in each bin for a 
uniform configuration (A^/2 = 245). The largest deviation 
from the average is 40 > y/N « 22, i.e. twice the statis- 
tic fiuctuation for a random distribution of particles. Main 
panel: The number of pairs of particles at a relative distance 
T'i < r < ri+i for the MD configuration in the inset, with 
To = 0, ri/a — 1 and r^+i — = a/10 for i > 1. The his- 
togram shows a large maximum corresponding to the attrac- 
tive range 2 < r/a < 2.2, a broad maximum around r/a = 1.2 
and a small number of pairs at the hard-core distance r = a. 
Therefore, the preferred relative distance for pairs of particles 
within the soft-core range 1 < r/a < 2 is, for this configura- 
tion, r/a ~ 1.2. 



To limit the computational effort, we consider 9 x 10^ 
wavevectors with modulus q < 100/a, that is much larger 
than the wavevectors of the largest peak of the crys- 
tal structure factor, q « 27r/(a/2), corresponding to the 
hard-core radius (Fig. |ll|c). Three different regimes can 
be distinguished in the example in Fig. |ll|. 

(i) A short-time regime "A", in which S{q,t) ^ 0(1) 
for any q (q = 1/a and q « 12/a w An/a are shown 



in Fig. 
in Fig. 



lla,b). Averaged on this interval (curve A 
Uc) the S{q) is fluid-like. 



(ii) An intermediate-time regime "B", in which S{q,t) 
for q = 1/a has an increase, but has no increase for 
q « An /a. Averaged on this interval (curve B in 
Fig. p]c) the S{q) is fluid-like, but with an increase 
for g — > 0. This increase indicates an increase of 
Kt, according to the equation 



ksTpRT — lim S{q) 



(13) 



where we use the Eq. (|^) and the deflnition S{q) = 
1 + ph{q). The increase of Kt is associated with 
the phase separation into two fluids with different 
densities. 

To help visualize the phase separation, in Fig. |l^ 
we show the three planar projections of the 3D con- 
figuration corresponding to the largest peak in the 
time interval "B" for g = 1/a (Fig. pT|a) . By divid- 
ing the box into two equal parts, the histograms of 
the number of particles in each part (Fig. |l^) show 
a separation in density approximately at half box- 
length, corresponding to g = An{p/Ny/^ « 1/a for 
p = 0.27/a^ and N = 490, in agreement with the 
peak at q — 1/a for the curve B in Fig. [Tl|c. In 
each projection, it is possible to see regions of high 
density and low density (Fig. p^ ). 

To better quantify the phase separation occurring 
in the configuration in Fig. |l^, we present (main 
panel in Fig. n3) the histogram of the number of 
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FIG. 13: Spatial distribution of pairs of particles at various 
distances for the MD configuration shown in Fig. The 
radius of the particles is not in scale with the distances. In 

(a) the darker particles are farther away from the observation 
point. Bonds connect particles at the hard-core distance 1 < 
r/a < 1.1 in panel (a); at distance 1.1 < r/a < 1.2 in panel 

(b) ; at distance 1.3 < r/a < 1.4 in panel (c); at distance 
1.5 < r/a < 1.6 in panel (d). The non-uniform distribution 
of bonds is clearly seen in panel (b) . 




FIG. 14: The last MD configuration in the time series in 
Fig. |l^. A crystal nucleus surrounded by gas is clearly seen. 
Bonds connect particles at distance r/a < 1.1. The 3D per- 
spective is given as in Fig. ^ The radius of particles is not in 
scale with the distances. 



the wavevectors corresponding to the peaks in the crystal 
S{q) and by a large step-like decrease of energy. 



C. The phase diagram and the finite size effect 



pairs of particles at a relative distance < r < 
Ti+i, where r^+i — rt — a/10. The histogram has 
a broad maximum around the distance r/a — 1.2 
in the soft-core range, showing that there exists 
a subset of pairs of particles that are at a pre- 
ferred distance 1.1 < r/a < 1.2. This subset is 
the HDL that has a non-uniform distribution over 
space (Fig. [T^b), consistent with the phase separa- 
tion. 

(iii) A long-time regime "C" , in which is S{q, t) ~ 0{N) 
for q « An /a, revealing the crystal nucleation pro- 
cess. The S{q) averaged over this time interval 
(curve C in Fig. pT|c) is solid-like. In the same 
interval, S{q,t) for q ~ \/a has a large increase, 
corresponding to the increase of Kt {S{q) increases 
for (7 — !■ 0), which is consistent with the phase sep- 
aration between the fluid and the crystal. As an 
example, in Fig. |lj we show the last configuration 
of the time series in Fig. |ri|, where the crystal struc- 
ture, already observed in Fig. ^, is clearly seen. 

The example in Fig. [Tl| b shows the formation of a high 
density fluid phase within the time interval "B" , followed 
by the nucleation of the crystal phase. The onset of the 
nucleation is marked by a large increase of S{q,t) for all 



By repeating the analysis described above for all the 
simulations inside the region with nucleation — and dis- 
carding the data corresponding to the formation of the 
nucleus — it is possible to calculate the state points cor- 
responding to the metastable fluid phase. The phase di- 
agram in Fig. |l^ is based on averages over a total of 
10^-10^ configurations in the fluid phase, accumulated 
in independent runs. 

For completeness we recall here the main features of 
the phase diagram in Fig. |l^ and presented in Ref. [ pT| . 
The (mechanically unstable) region at high p for T < 
0.67J7yi//cB, where P decreases for increasing p, denotes 
the coexistence of gas and HDL. The unstable region at 
low p for T < 0.603UA/kB (inset in Fig. ^ denotes 
the coexistence of gas and LDL. The coexistence lines 
are obtained by using the Maxwell construction of the 
equal areas 0], suggesting the presence of a gas-LDL- 
HDL triple point. 

By definition, the spinodal lines (limit of stability of 
each phase with respect to the coexisting phase) meet 
the coexistence lines in a critical point. Therefore, by in- 
terpolation we estimate the gas-LDL critical point Ci at 
keTi/UA = 0.603±0.003, a^i = O.lOiO.Ol, a^Pi/UA = 
0.0171 ± 0.0005 and, the gas-HDL critical point C2 
at kBT2/UA = 0.665 ± 0.005, a^p2 = 0.306 ± 0.020, 
a^P2/UA = O.lOibO.Ol. These values are consistent with 




FIG. 15: Comparison between MD simulations for — 490 
(full diamonds) and iV = 1728 (circles) at ksT/UA = 0.595. 
The results for the two sizes are very close. For comparison, 
we include also the calculations for N = 490 at ksT /Ua = 
0.60 and 0.59 (upper open diamonds and lower open dia- 
monds, respectively) and the interpolation at ksT/UA = 
0.595 (dashed line) between these two isotherms, showing the 
presence of two regions with negatively sloped isotherms. The 
points calculated for A'' = 1728 are also consistent with this 
interpolation, suggesting that the finite size effect between 
A'' = 490 and A'' = 1728 is small. Errors, where not shown, 
are smaller than the symbols. 



the linear interpolations of the MD isotherms (Fig. |T^. 

The phase diagram resulting from the MD calculations 
is, as expected, in agreement with the time-dependent 
analysis of the structure factor presented above. For 
example, the case presented in Fig. 11 corresponds to 
a state point inside the gas-HDL coexistence region at 
a density higher then the crystal nucleation density for 
T = 0.62UA/kB- The nucleation of the (metastable) 
HDL phase is thus followed by the crystal nucleation. 

To estimate the finite size effect in our calculations, 
we compare the results for N — 490 and N = 1728 for 
an isotherm below both critical points (Fig. |l^). The 
calculations do not show any relevant finite size effect, 
suggesting that the MD results for = 490 are reliable. 



V. THE RADIAL DISTRIBUTION FUNCTION 
ANALYSIS 

The interpretation of the HNC instability line is quali- 
tatively consistent with the MD spinodal line for the cor- 
responding set of the potential's parameters. The pro- 
jection of the MD spinodal line in the T-p plane (not 
shown) has the same characteristics of the HNC instabil- 
ity line, with two local maxima and one local minimum. 
In both approaches, the high-p local maximum occurs at 
a temperature higher then the temperature of the low-p 
maximum and the presence of a triple point is suggested 



FIG. 16: Comparison between the g{r) calculated in the HNC 
approximation (dotted line) and by MD simulations (solid 
line). As an example, we present the calculations at T = 
0.6iUA/kB for densities (from top to bottom) a^p = 0.066, 
0.154, 0.267. For clarity a constant value is added to the first 
two curves (6.5 and 3, respectively). The two independent 
calculations are very close at intermediate densities. At large 
r all the curves go to 1. 



by the presence of the local minimum. 

The quantitative HNC predictions for the locations of 
the two critical points are, as expected, only partially 
consistent with the MD results. It is remarkable that the 
HNC estimates of the density of the low-p local maxi- 
mum (p w and the temperature of the high-p local 
maximum (T « 0.65J7a/^_b) of the instability line are 
close to the corresponding MD results for Ci and C2, 
respectively. 

An estimate of the agreement between the two meth- 
ods can be evaluated by comparing the calculations for 
g{r) within the two approaches (Fig. |l^). In contrast 
with what could be suggested by the nature of the HNC 
approximation — i.e. the under-estimate of the indirect 
correlation — the agreement is better at intermediate p 
than at low p (Fig. |l^). In particular, at low p the 
HNC approximation underestimates the probability of a 
particle penetrating the soft-core or entering the attrac- 
tive well. At higher p, instead, the estimates of the g{r) 
within the two approaches are almost indistinguishable. 

The g{r) of the low-p fluid is characterized by a large 
peak at r = 6 corresponding to the shortest attractive 
distance. As a consequence of the increase of the den- 
sity, the peak at the hard-core distance r — a increases 
while the peak at r = 6 decreases, and additional peaks 
at r/a = 3, 4, ... appear. In Fig. |l^ we present the 
calculation of the g{r) for the gas phase, the gas-HDL 
coexisting region and the HDL phase. In particular, by 
combining the radial distribution functions evaluated in 
each pure phase, we can estimate the composition of the 
mixed phase. For example, at T = 0.64Ua/ ks the ra- 
dial distribution function calculated at po — 0.302/a'^ is 
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FIG. 17: The radial distribution function g{r) calculated from 
the MD results at T = 0.64(7A/fcB for densities a^p = 0.066, 
0.100, 0.132, 0.154, 0.186, 0.223, 0.267, 0.302, 0.322, 0.333, 
0.349, 0.358. With increasing p, the peak at r = a increases 
and the peak at r = b decreases, while more peaks appear at 
larger r. 
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FIG. 19: Inset: the MD resuhs at T = 0.64;7a/A:b for the cu- 
mulative number of particles A7V within the repulsive range 
r < b (circles) and within the attractive range b < r < c 
(squares), increasing linearly with p in this range of densities; 
the linear fit of the data gives a slope 20.2 ± 0.5 for the solid 
line and a slope 22.1 ± 0.6 for the dashed line. Main panel: 
the data in the inset plotted one versus the other to show that 
within the mechanically unstable region (p > 0.267/a^) AiV 
increases approximately in the same way with p both within 
the attractive and within the repulsive range (the dashed line 
is a linear fit of all the data with slope 1.08 ± 0.06) and in- 
creases faster within the attractive range at lower densities 
(the solid line is a quadratic fit). 



r (Fig. |l9|). This analysis reveals that the number of 
particles AiV within the repulsive range and within the 
attractive range increases linearly with p (inset Fig. [l9| ) 
and that the increase is faster within the attractive range 
(Fig. |l^), for the densities we studied. In particular, the 
number of particles within the attractive range b < r < c 
increases from 2.5 to 9 < Em=i<ryiO = 19.2 esti- 
mated for the artificial crystal (Table O). 



FIG. 18: The radial distribution function calculated from the 
MD simulations at T = 0.64UA/kB and p = 0.302/a^ (open 
circles) is compared with the composition Xigi{r) + X2g2{r) 
(solid line), where gi{r) is the radial distribution function for 
the pure gas phase (at p — 0.223/a'^) and g2{r) is for the pure 
HDL phase (at p = 0.349/a^) at the same T, with Xi = 0.3 
and X2 = 0.7. 



go{r) ~ Xigi{r) -f X2g2{r), where gi{r) and 32 (r) are 
the radial distribution functions at the same T and at 
Pi — 0.223/a'^ and p2 = 0.349/a^, respectively, with 
Xi = 0.3 and X2 = 0.7 (Fig. |l^, revealing that the 
system is composed approximately by 30% of gas and 
70% of HDL. 

From the Eq. (|^), by using the g{r) calculated from 
the MD simulations, we evaluate the average number 
of particles N{r) = J dN{r) within a sphere of radius 



VI. ABSENCE OF A DENSITY ANOMALY 

In Ref. |jll| it has been noted that the possibility of 
a second fluid-fluid critical point is not necessarily re- 
stricted to systems with a density anomaly, at least from 
a theoretical point of view. Here we present the explicit 
thermodynamic calculations for this result. 

The defining relation for the density anomaly is given 

by 
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df 
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FIG. 20: The potential energy density U/N, calculated by 
MD simulations, as a function of the density p for tempera- 
tures (bottom to top) ksT/UA = 0.60, 0.61, 0.63, 0.64, 0.66, 
0.67, 0.70. The symbols represent the MD calculations, with 
errors smaller than the symbol's size. The lines represent the 
cubic fit of the calculations with the parameters in Table III 



for the Maxwell relation 

dV 

df 

where S is the entropy. Since 

dV 
dP 



dS_ 



< 



(16) 



(17) 



holds for a mechanically stable phase, Eq. (|lj) can be 
rewritten as 



dS_ 

dV 



< 



(18) 



From the differential expression of the thermodynamic 
potential at constant T, we know that 



TdS = dE + PdV , 



(19) 



where E = U + K is the total energy, with U and K total 
potential and kinetic energy, respectively. Therefore, it 
is 
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T dV 
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(20) 



at constant T and we can rewrite the density anomaly 
condition in Eq. (|l^) as 



dV 



P{V, T) < 



(21) 



at constant T. 

To calculate the left-hand side of Eq. (pT]), we need 
to evaluate {dU/dV)T- In Fig. 20, we show our MD 



TABLE III: Parameters for the cubic fit U/N = ao + aip + 
a2p^ + aap^ of the MD calculations for the potential energy 
density U/N in Fig. ^ for different temperatures. The errors 
on the fitting parameters are on the last decimal digit. 
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FIG. 21: The derivative (dU/dV)T calculated by using the 
cubic expression in Fig. bol with the parameters in Table [II , 
as function of the specific volume V/N for temperatures (top 
to bottom) ksT/UA = 0.60, 0.61, 0.63, 0.64, 0.66, 0.67, 0.70. 
The derivative, in the considered range of V/N, is always 
larger than 0.025f7A/a^ (bottom horizontal line). Inset: The 
same MD resuhs in Fig. M for U/N plotted as a function of 
V/N. The symbols are the same as in Fig. po| . 



calculation for U{p) at constant T. All the MD points 
can be fitted with a third degree polynomial in p. The 
fitting parameters are given in Table III and are used 
to calculate the derivative (dU / dV)T , shown in Fig. |2^. 
Our calculations show a potential energy U increasing 
with V (inset Fig. |2^), with a derivative always positive, 
thus wherever P is positive, the condition in Eq. ( ^ ) is 
not satisfied and there is no density anomaly. 

In the region where P < (at low T and small V), 
the derivative {dU / dV)T rapidly increases in such a way 
that Eq. (|l|) is never satisfied. Particularly, in the 
range of volumes considered, it is always {dU/dV)T — 
0.025f7A/a3 > 0, where P = ~0.025UA/a^ is the mini- 
mum pressure, reached for T — O.GC/^/fcs and V/N — 
3.31a'^ (Fig. |lO|). These results suggest that the density 
anomaly is ruled out for this choice of parameters. At 
this stage is not clear if it is ruled out for any choice of 
parameters for our potential in 3D (see Ref. [p^). 
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VII. SUMMARY AND CONCLUSIONS 

We analyzed the phase diagram of a soft-core poten- 
tial, similar to potentials used in systems such as pro- 
tein solutions, colloids, melts and in pure systems such 
as liquid metals. We use two independent numerical 
methods, integral equations in the HNC approximation 
and MD simulations. The comparison of the HNC re- 
sults with previously-proposed soft-core potentials sug- 
gests that the system has two fluid-fluid phase transi- 
tions for an appropriate choice of parameters — energy 
and width — of the repulsive soft-core. We select a set of 
potential parameters with a narrow attractive well that 
gives a HNC instability line with two maxima and sug- 
gests the presence of two critical points. 

The MD analysis shows, in agreement with the pre- 
vious results for potentials with a short range attraction 
p2| , a phase diagram with no stable liquid phase. We an- 
alyze the crystal structure, characterized by the competi- 
tion between the attractive interaction at distance r — b 
and the repulsive interaction at r = a < b. We show 
that the crystal, with 8-fold and 12-fold symmetries, is 
independent on the density, within the considered range 
of densities. 

Hence, we study the metastable liquid phase, at tem- 
peratures above and below the line of spontaneous crystal 
nucleation. We find two liquids in the supercooled phase, 
the LDL and the HDL, with two fluid-fluid transitions 
ending in two critical points, the gas-LDL critical point 
Ci and the gas-HDL critical point C2, as already shown 



in Ref. [|ll|. Here we improve our estimate of the phase 
diagram and we verify the absence of relevant finite size 
effects in the MD results. 

We compare these results with the HNC calculations, 
concluding that the HNC approximation underestimates 
the effect of the attractive interaction and overestimates 
the effect of the repulsive interaction at low p, and is in 
good agreement with the MD results at intermediate p. 

Finally, by explicit calculation, we show that the condi- 
tion for the density anomaly is never satisfied in the range 
of T and V considered here, as announced in Ref. [ pT| . 
Our results suggest that the density anomaly is always 
ruled out for this choice of potential parameters. 

In conclusion, the results of this paper evoke an intrigu- 
ing relation between the absence of the density anomaly 
and the presence of a single crystal phase, with higher 
density than the liquid phases, in systems with two 
fluid-fluid phase transitions. This relation, that deserves 
greater investigation, is consistent with the fact that the 
substances with the density anomaly, and a hypothesized 
second liquid- liquid critical point, have more than one 
crystal structure, as in the case of water or carbon or 
silica. 
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